rm(list=ls(all=TRUE))
graphics.off()
library(rdd)
library(xtable)
library(dplyr)

setwd() # set working directory
load("data/data_alt_boot.rdata")

# limit to parties with at least one elected 
data_alt_elected <-
  data_alt %>%
  group_by(year, muncpr, candpart) %>%
  summarise(represented = sum(elected)>0)

data_alt <- 
  left_join(data_alt, data_alt_elected) %>%
  filter(represented == 1)

#drop all from closed and split lists 

data_alt <- 
  data_alt %>%
  filter(openlist == 1 & splitlist == 0)

# Create clusters of year, municipality, and party

clusters           <- data.frame(unique(cbind(data_alt$year,data_alt$muncpr,data_alt$candpart)))
colnames(clusters) <- c("year", "muncpr", "candpart")
clusters$cluster   <- 1:nrow(clusters)
data_alt$year        <- as.factor(data_alt$year)
data_alt$muncpr      <- as.factor(data_alt$muncpr)
data_alt <- left_join(data_alt, clusters, by = c("year", "muncpr", "candpart"))

#drop if the marginal loser never wins in sim or the marginal winners always win
data_alt <- 
  data_alt %>%
  filter(marg.los != 0 & marg.win != 1)

#drop candidates that always or never win

data_alt <- 
  data_alt %>%
  filter(simscore > 0 & simscore < 1)

#if cases are extremely close they might have wrong signed simscore. Recode these to arbitrarily small simscore
data_alt$margsim[data_alt$margsim <= 0 & data_alt$elected == 1] <-  0.0001
data_alt$margsim[data_alt$margsim >= 0 & data_alt$elected == 0] <- -0.0001

#create variables for belonging to main left wing parties (including Social Liberals)
#and personal votes as proportion of list votes and municipality total votes
data_alt$mainleft      <- data_alt$candpart %in% c("A", "F", "B", "Oe")
data_alt$share_of_party <- data_alt$candvote/data_alt$ptotvote
data_alt$share_of_mun  <- data_alt$candvote/data_alt$mtotvavo

data_2005 <- 
  data_alt %>%
  filter(year == 2005)

data_2009 <- 
  data_alt %>%
  filter(year == 2009)

# Estimate and plot RD for having incumbency in 2009
# Find Imbens-Kalyanaraman bandwidth
IK.bw_2005  <- IKbandwidth(data_2005$margsim, data_2005$electedt1)
IK.bw2_2005 <- IKbandwidth(data_2005$margsim, data_2005$rerun)

##rerunning and elected 
rdest_2005    <- RDestimate(electedt1 ~ margsim, data = data_2005)
rdest.cl_2005 <- RDestimate(electedt1 ~ margsim, data = data_2005, cluster = data_2005$cluster)

##rerunning  
rdest1_2005    <- RDestimate(rerun ~ margsim, data = data_2005)
rdest1.cl_2005 <- RDestimate(rerun ~ margsim, data = data_2005, cluster = data_2005$cluster)

##Summarize jumps in table

out1_2005 <- cbind(rdest_2005$est,
                   rdest_2005$se,
                   rdest.cl_2005$se,
                   rdest_2005$bw,
                   rdest_2005$obs)
colnames(out1_2005) <- c("Estimate", "Std. Error", "Rob. Std. Error", "Bandwidth", "Observations")

out2_2005 <- cbind(rdest1_2005$est,
              rdest1_2005$se,
              rdest1.cl_2005$se,
              rdest1_2005$bw,
              rdest1_2005$obs)
colnames(out2_2005) <- c("Estimate", "Std. Error", "Rob. Std. Error", "Bandwidth", "Observations")

digits_tab <- rbind(matrix(3, nrow = 4, ncol = 7), rep(0,7))


# 2013 --------------------------------------------------------------------

# Estimate and plot RD for having incumbency in 2013
# Find Imbens-Kalyanaraman bandwidth
IK.bw_2009  <- IKbandwidth(data_2009$margsim, data_2009$electedt1)
IK.bw2_2009 <- IKbandwidth(data_2009$margsim, data_2009$rerun)

##rerunning and elected 
rdest_2009    <- RDestimate(electedt1 ~ margsim, data = data_2009)
rdest.cl_2009 <- RDestimate(electedt1 ~ margsim, data = data_2009, cluster = data_2009$cluster)

##rerunning  
rdest1_2009    <- RDestimate(rerun ~ margsim, data = data_2009)
rdest1.cl_2009 <- RDestimate(rerun ~ margsim, data = data_2009, cluster = data_2009$cluster)

##Summarize jumps in table

out1_2009 <- cbind(rdest_2009$est,
                   rdest_2009$se,
                   rdest.cl_2009$se,
                   rdest_2009$bw,
                   rdest_2009$obs)
colnames(out1_2009) <- c("Estimate", "Std. Error", "Rob. Std. Error", "Bandwidth", "Observations")

out2_2009 <- cbind(rdest1_2009$est,
                   rdest1_2009$se,
                   rdest1.cl_2009$se,
                   rdest1_2009$bw,
                   rdest1_2009$obs)
colnames(out2_2009) <- c("Estimate", "Std. Error", "Rob. Std. Error", "Bandwidth", "Observations")

digits_tab <- rbind(matrix(3, nrow = 4, ncol = 7), rep(0,7))


# Print tables
xtable(cbind(t(out1_2005),t(out2_2005)), digits = digits_tab)
xtable(cbind(t(out1_2009),t(out2_2009)), digits = digits_tab)

graphics.off()
quartz(w = 6, h = 5)
par(mar=c(5,7,3,1))
# create plot of effects over year
plot(x = c(rdest_2005$est[1], rdest1_2005$est[1]),
     y = 2:1 + 0.1,
     axes = FALSE,
     ylim = c(0.5, 2.5),
     xlim = c(-.1, 0.35),
     xlab = "Effect size",
     ylab = "",
     pch  = 19)

segments(x0 = c(rdest_2005$est[1] + rdest_2005$se[1]*qnorm(0.025),
                rdest1_2005$est[1] + rdest1_2005$se[1]*qnorm(0.025)),
         y0 = 2:1 + 0.1,
         x1 = c(rdest_2005$est[1] + rdest_2005$se[1]*qnorm(0.975),
                rdest1_2005$est[1] + rdest1_2005$se[1]*qnorm(0.975)),
         y1 = 2:1 + 0.1,
         lwd = 3)

# Add 2009
points(x = c(rdest_2009$est[1], rdest1_2009$est[1]),
     y = 2:1 - 0.1,
     col  = "grey50",
     pch  = 19)

segments(x0 = c(rdest_2009$est[1] + rdest_2009$se[1]*qnorm(0.025),
                rdest1_2009$est[1] + rdest1_2009$se[1]*qnorm(0.025)),
         y0 = 2:1 - 0.1,
         x1 = c(rdest_2009$est[1] + rdest_2009$se[1]*qnorm(0.975),
                rdest1_2009$est[1] + rdest1_2009$se[1]*qnorm(0.975)),
         y1 = 2:1 - 0.1,
         lwd = 3,
         col = "grey50")

legend("topright", 
       pch = c(19, 19), 
       c("Advantage in 2009",
         "Advantage in 2013"),
       col = c("black", "grey"),
       bty = "n")
axis(1, at = seq(-0.1, 0.35, 0.05))
par(las = 1)

axis(2, at = 2:1, 
     labels = c("p(rerunning \nand elected)",
                "p(rerunning)"),
     col = "white")



